#pragma once
#include <fstream>
#include <cmath>
#include <complex>
#include <iostream>
#include <iomanip>
#include <limits>
#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <fcntl.h>
#include <string.h>
#include <ctype.h>

namespace mylib
{
  template <class T>
  inline T SQR(const T a) { return a * a; }

  template <class T>
  inline const T &MAX(const T &a, const T &b)
  {
    return b > a ? (b) : (a);
  }

  inline float MAX(const double &a, const float &b)
  {
    return b > a ? (b) : float(a);
  }

  inline float MAX(const float &a, const double &b)
  {
    return b > a ? float(b) : (a);
  }

  template <class T>
  inline const T &MIN(const T &a, const T &b)
  {
    return b < a ? (b) : (a);
  }

  inline float MIN(const double &a, const float &b)
  {
    return b < a ? (b) : float(a);
  }

  inline float MIN(const float &a, const double &b)
  {
    return b < a ? float(b) : (a);
  }

  template <class T>
  inline T SIGN(const T &a, const T &b)
  {
    return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a);
  }

  inline float SIGN(const float &a, const double &b)
  {
    return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a);
  }

  inline float SIGN(const double &a, const float &b)
  {
    return (float)(b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a));
  }

  template <class T>
  inline void SWAP(T &a, T &b)
  {
    T dum = a;
    a = b;
    b = dum;
  }

  template <class T>
  class vector
  {
  private:
    int n; //矢量长度
    T *v;  //数据
  public:
    vector() : n(0), v(NULL) //空的初始化
    {
    }
    explicit vector(int n) : n(n), v(n > 0 ? new T[n] : NULL) //空的初始化,禁止隐式转换
    {
    }
    vector(int n, const T &a) : n(n), v(n > 0 ? new T[n] : NULL) //用常数来初始化
    {
      for (int i = 0; i < n; i++)
        v[i] = a;
    }
    vector(int n, const T *a) : n(n), v(n > 0 ? new T[n] : NULL) //用一维数组来初始化
    {
      for (int i = 0; i < n; i++)
        v[i] = *a++;
    }
    vector(const vector &rhs) : n(rhs.n), v(n > 0 ? new T[n] : NULL) //用该类来初始化
    {
      for (int i = 0; i < n; i++)
        v[i] = rhs[i];
    }

    vector<T> &operator=(const vector<T> &rhs) // 重载赋值操作
    {
      if (this != &rhs)
      {
        if (n != rhs.n)
        {
          if (v != NULL)
            delete[](v);
          n = rhs.n;
          v = n > 0 ? new T[n] : NULL;
        }
        for (int i = 0; i < n; i++)
          v[i] = rhs[i];
      }
      return *this;
    }
    typedef T value_type; //

    inline T &operator[](const int i) //重载[]算符
    {
      return v[i];
    }

    inline const T &operator[](const int i) const
    {
      return v[i];
    }

    inline int size() const
    {
      return n;
    }

    void resize(int newn)
    {
      if (newn != n)
      {
        if (v != NULL)
          delete[](v);
        n = newn;
        v = n > 0 ? new T[n] : NULL;
      }
    }
    void assign(int newn, const T &a)
    {
      if (newn != n)
      {
        if (v != NULL)
          delete[](v);
        n = newn;
        v = n > 0 ? new T[n] : NULL;
      }
      for (int i = 0; i < n; i++)
        v[i] = a;
    }

    void show(int nn)
    {
      int newn = (nn > n ? n : nn);
      if (v != NULL)
      {
        //          std::cout<<"v[i] = [";
        for (int i = 0; i < newn; i++)
        {
          std::cout << std::setw(5) << v[i] << ","; // c++ iostream
        }
        //          std::cout<<"]";
        std::cout << std::endl;
      }
    }

    vector<T> operator*(const vector<T> &a)
    {
      int newn = std::min(a.n, n);
      vector<T> c(newn);

      for (int i = 0; i < newn; i++)
        c[i] = v[i] * a[i];

      return c;
    }

    vector<T> operator*(const T &a)
    {
      vector<T> c(n);

      for (int i = 0; i < n; i++)
        c[i] = v[i] * a;

      return c;
    }

    vector<T> operator+(const vector<T> &a)
    {
      int newn = std::min(a.n, n);
      vector<T> c(newn);

      for (int i = 0; i < newn; i++)
        c[i] = v[i] + a[i];

      return c;
    }

    vector<T> operator-(const vector<T> &a)
    {
      int newn = std::min(a.n, n);
      vector<T> c(newn);

      for (int i = 0; i < newn; i++)
        c[i] = v[i] - a[i];

      return c;
    }

    T dot(const vector<T> &a)
    {
      T sum = T(0);
      int newn = std::min(a.n, n);

      for (int i = 0; i < newn; i++)
      {
        sum = sum + a[i] * v[i];
      }

      return sum;
    }

    T normal()
    {
      T sum = T(0);
      for (int i = 0; i < n; i++)
      {
        sum = sum + v[i] * v[i];
      }
      return sqrt(sum);
    }

    ~vector()
    {
      if (v != NULL)
        delete[](v);
    }

  }; // class vector

  template <class T>
  class matrix
  {
  private:
    int n;
    int m;
    T **v;

  public:
    matrix() : n(0), m(0), v(NULL) {}
    matrix(int n, int m) : n(n), m(m), v(n > 0 ? new T *[n] : NULL) //零矩阵
    {
      int i, nel = m * n;
      if (v)
        v[0] = nel > 0 ? new T[nel] : NULL;
      for (i = 1; i < n; i++)
        v[i] = v[i - 1] + m;
    }

    matrix(int n, int m, const T &a) : n(n), m(m), v(n > 0 ? new T *[n] : NULL) // 用常数初始化
    {
      int i, j, nel = m * n;
      if (v)
        v[0] = nel > 0 ? new T[nel] : NULL;
      for (i = 1; i < n; i++)
        v[i] = v[i - 1] + m;
      for (i = 0; i < n; i++)
        for (j = 0; j < m; j++)
          v[i][j] = a;
    }

    matrix(int n, int m, const T *a) : n(n), m(m), v(n > 0 ? new T *[n] : NULL) // 用数组初始化
    {
      int i, j, nel = m * n;
      if (v)
        v[0] = nel > 0 ? new T[nel] : NULL;
      for (i = 1; i < n; i++)
        v[i] = v[i - 1] + m;
      for (i = 0; i < n; i++)
        for (j = 0; j < m; j++)
          v[i][j] = *a++;
    }

    matrix(const matrix &rhs) : n(rhs.n), m(rhs.m), v(n > 0 ? new T *[n] : NULL) // 复制构造函数
    {
      int i, j, nel = m * n;
      if (v)
        v[0] = nel > 0 ? new T[nel] : NULL;
      for (i = 1; i < n; i++)
        v[i] = v[i - 1] + m;
      for (i = 0; i < n; i++)
        for (j = 0; j < m; j++)
          v[i][j] = rhs[i][j];
    }

    matrix<T> &operator=(const matrix<T> &rhs) // 重载赋值运算符
    {
      if (this != &rhs)
      {
        int i, j, nel;
        if (n != rhs.n || m != rhs.m)
        {
          if (v != NULL)
          {
            delete[](v[0]);
            delete[](v);
          }
          n = rhs.n;
          m = rhs.m;
          v = n > 0 ? new T *[n] : NULL;
          nel = m * n;
          if (v)
            v[0] = nel > 0 ? new T[nel] : NULL;
          for (i = 1; i < n; i++)
            v[i] = v[i - 1] + m;
        }
        for (i = 0; i < n; i++)
          for (j = 0; j < m; j++)
            v[i][j] = rhs[i][j];
      }
      return *this;
    }

    typedef T value_type; //

    inline T *operator[](const int i) // 下标第i行
    {
      return v[i];
    }

    inline const T *operator[](const int i) const // 下标第i行
    {
      return v[i];
    }

    inline int nrows() const
    {
      return n;
    }

    inline int ncols() const
    {
      return m;
    }

    vector<T> selcol(const int i,const int i_start, const int i_end)
    {
      vector<T> x(i_end-i_start);
      for (int j = 0; j < x.size(); j++)
      {
        x[j] = v[j+i_start][i];
      }
      return x;
    }
    void setcol(const int i,const int i_start,vector<T> x)
    {
      for (int j = 0; j < x.size(); j++)
      {
        v[j+i_start][i]=x[j] ; 
      }
    }

    vector<T> selrow(const int i,const int i_start, const int i_end)
    {
      vector<T> x(i_end-i_start);
      for (int j = 0; j < x.size(); j++)
      {
        x[j] = v[i][j+i_start];
      }
      return x;
    }

    void setrow(const int i,const int i_start,vector<T> x)
    {
      for (int j = 0; j < x.size(); j++)
      {
        v[i][i_start+j]=x[j] ; 
      }
    }

    matrix<T> selblock(const int i, const int j)
    {
      int m = j - i;
      matrix<T> blockA(m, m);

      for (int ii = i; ii < j; ii++)
      {
        for (int jj = i; jj < j; jj++)
        {
          blockA[ii, jj] = v[ii][jj];
        }
      }
      return blockA;
    }

    void resize(int newn, int newm) // 改变大小，内容会丢失
    {
      int i, nel;
      if (newn != n || newm != m)
      {
        if (v != NULL)
        {
          delete[](v[0]);
          delete[](v);
        }
        n = newn;
        m = newm;
        v = n > 0 ? new T *[n] : NULL;
        nel = m * n;
        if (v)
          v[0] = nel > 0 ? new T[nel] : NULL;
        for (i = 1; i < n; i++)
          v[i] = v[i - 1] + m;
      }
    }
    void assign(int newn, int newm, const T &a) // 改变大小，用常数赋值
    {
      int i, j, nel;
      if (newn != n || newm != m)
      {
        if (v != NULL)
        {
          delete[](v[0]);
          delete[](v);
        }
        n = newn;
        m = newm;
        v = n > 0 ? new T *[n] : NULL;
        nel = m * n;
        if (v)
          v[0] = nel > 0 ? new T[nel] : NULL;
        for (i = 1; i < n; i++)
          v[i] = v[i - 1] + m;
      }
      for (i = 0; i < n; i++)
        for (j = 0; j < m; j++)
          v[i][j] = a;
    }

    matrix<T> dot(matrix<T> &a)
    {
      matrix<T> sum(n, a.m);
      if (m != a.n)
      {
        std::cout << "a.row ne b.col" << std::endl;
        return sum;
      }

      for (int i = 0; i < n; i++)
        for (int j = 0; j < a.m; j++)
          for (int k = 0; k < m; k++)
          {
            sum[i][j] += v[i][k] * a[k][j];
          }

      return sum;
    }

    void show(int nn, int mm)
    {
      int newn = (nn > n ? n : nn);
      int newm = (mm > m ? m : mm);

      std::cout << "v[i][j]=" << std::endl;
      for (int i = 0; i < newn; i++)
      {
        for (int j = 0; j < newm; j++)
        {
          std::cout << std::setw(5) << v[i][j] << ",";
        }
        std::cout << std::endl;
      }
    }

    ~matrix()
    {
      if (v != NULL)
      {
        delete[](v[0]);
        delete[](v);
      }
    }
  };

  typedef int Int; // 32 bit integer
  typedef unsigned int Uint;

#ifdef _MSC_VER
  typedef __int64 Llong; // 64 bit integer
  typedef unsigned __int64 Ullong;
#else
  typedef long long int Llong; // 64 bit integer
  typedef unsigned long long int Ullong;
#endif

  typedef char Char; // 8 bit integer
  typedef unsigned char Uchar;

  typedef double Doub; // default floating type
  typedef long double Ldoub;

  typedef std::complex<double> Complex; // default complex type

  typedef bool Bool;

  // NaN: uncomment one of the following 3 methods of defining a global NaN
  // you can test by verifying that (NaN != NaN) is true

  // static const Doub NaN = numeric_limits<Doub>::quiet_NaN();

  // Uint proto_nan[2]={0xffffffff, 0x7fffffff};
  // double NaN = *( double* )proto_nan;

  // Doub NaN = sqrt(-1.);

  // vector types

  typedef const vector<Int> VecInt_I;
  typedef vector<Int> VecInt, VecInt_O, VecInt_IO;

  typedef const vector<Uint> VecUint_I;
  typedef vector<Uint> VecUint, VecUint_O, VecUint_IO;

  typedef const vector<Llong> VecLlong_I;
  typedef vector<Llong> VecLlong, VecLlong_O, VecLlong_IO;

  typedef const vector<Ullong> VecUllong_I;
  typedef vector<Ullong> VecUllong, VecUllong_O, VecUllong_IO;

  typedef const vector<Char> VecChar_I;
  typedef vector<Char> VecChar, VecChar_O, VecChar_IO;

  typedef const vector<Char *> VecCharp_I;
  typedef vector<Char *> VecCharp, VecCharp_O, VecCharp_IO;

  typedef const vector<Uchar> VecUchar_I;
  typedef vector<Uchar> VecUchar, VecUchar_O, VecUchar_IO;

  typedef const vector<Doub> VecDoub_I;
  typedef vector<Doub> VecDoub, VecDoub_O, VecDoub_IO;

  typedef const vector<Doub *> VecDoubp_I;
  typedef vector<Doub *> VecDoubp, VecDoubp_O, VecDoubp_IO;

  typedef const vector<Complex> VecComplex_I;
  typedef vector<Complex> VecComplex, VecComplex_O, VecComplex_IO;

  typedef const vector<Bool> VecBool_I;
  typedef vector<Bool> VecBool, VecBool_O, VecBool_IO;

  // matrix types

  typedef const matrix<Int> MatInt_I;
  typedef matrix<Int> MatInt, MatInt_O, MatInt_IO;

  typedef const matrix<Uint> MatUint_I;
  typedef matrix<Uint> MatUint, MatUint_O, MatUint_IO;

  typedef const matrix<Llong> MatLlong_I;
  typedef matrix<Llong> MatLlong, MatLlong_O, MatLlong_IO;

  typedef const matrix<Ullong> MatUllong_I;
  typedef matrix<Ullong> MatUllong, MatUllong_O, MatUllong_IO;

  typedef const matrix<Char> MatChar_I;
  typedef matrix<Char> MatChar, MatChar_O, MatChar_IO;

  typedef const matrix<Uchar> MatUchar_I;
  typedef matrix<Uchar> MatUchar, MatUchar_O, MatUchar_IO;

  typedef const matrix<Doub> MatDoub_I;
  typedef matrix<Doub> MatDoub, MatDoub_O, MatDoub_IO;

  typedef const matrix<Bool> MatBool_I;
  typedef matrix<Bool> MatBool, MatBool_O, MatBool_IO;

}